This notebook implements a GHSL settlement model algorithm as defined by the stage I of the Degree of Urbanisation (European Commission & Statistical Office of the European Union, 2021) and recommended by the UN STAT COM.
Model uses the population and built-up surface grid - European Commission, Joint Research Centre (JRC)
The method here implemented is described in detail bellow.
Delineate and classify settlement by typologies
GHSL data packages for settlement model:
# local population density (HDC: 1500 people/km2, MDC: 300 people/km2)
HDC_Pdens=1500
MDC_Pdens =300
# built-up area share > then
HDC_Bdens=500000 # (share of built-up area/km2: 0.5 => 500000 m2/km2)
# cluster populatin > then (HDC: 50000 people, MDC: 5000 people)
HDC_Pmin=50000
MDC_Pmin=5000
# connectivity clusters
# Tcon=1 for 4-connectivity clusters, Tcon=2 for 8-connectivity clusters
HDC_Tcon=1
MDC_Tcon=2
# gap filling (<15km2)
HDC_GAPfill=15
from xcube.core.store import new_data_store
import os
from glob import glob
from pathlib import Path
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import shape, box
import xarray as xr
import rioxarray as rxr
import rasterio
import rasterio.features
from rasterio.features import shapes
from skimage import measure, morphology
from skimage.morphology import square
import scipy
import ipywidgets as widgets
import matplotlib.pyplot as plt
import folium
from folium import plugins
/opt/conda/envs/edc-default-2022.10-14/lib/python3.9/site-packages/geopandas/_compat.py:112: UserWarning: The Shapely GEOS version (3.11.0-CAPI-1.17.0) is incompatible with the GEOS version PyGEOS was compiled with (3.10.3-CAPI-1.16.1). Conversions between both will be slow. warnings.warn(
Before starts, the reference year and the area of interest needs to be defined.
# Select reference year:
# from 1975 to 2020, 5 years interval:
# valid values are: 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020
REFyear=2020
# Select AOI by lat/lon coordinates - example for area of Poland, Czechia and Slovakia:
x1 = 10 # degree (E-W)
y1 = 47 # degree (N-S)
x2 = 25 # degree (E-W)
y2 = 55 # degree (N-S)
aoi_centre = [(y2+y1)/2,(x2+x1)/2]
print(aoi_centre)
bbox = x1, y1, x2, y2
print(bbox)
[51.0, 17.5] (10, 47, 25, 55)
geom = box(*bbox)
AOI = gpd.GeoDataFrame({"id":1,"geometry":[geom]})
AOI.set_crs(epsg=4326, inplace=True) # WGS84
AOI54009 = AOI.to_crs(crs='ESRI:54009')
from branca.element import Figure
fig = Figure(width=600, height=600)
m = folium.Map(aoi_centre, zoom_start=5, tiles='cartodbpositron') # ,width=500, height=500
lat_interval = 1
lon_interval = 1
# parallels:
for lat in range(-90, 91, lat_interval):
folium.PolyLine([[lat, -180],[lat, 180]], weight=0.5).add_to(m)
# meridianss:
for lon in range(-180, 181, lon_interval):
folium.PolyLine([[-90, lon],[90, lon]], weight=0.5).add_to(m)
AOI54009.explore(m=m,style_kwds=dict(color='blue',fillColor='None', opacity=0.7))
fig.add_child(m)
m
BTU_data_store = new_data_store("s3", root="xcube-dcfs/GHSL/GHS-BUILT-S/",
storage_options=dict(anon=True) # anon=True - anonymous access (default False)
)
BTU_data_store_info = list(BTU_data_store.get_data_ids())
BTU_data_store_info
['GHS_BUILT_S_E1975_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E1980_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E1985_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E1990_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E1995_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E2000_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E2005_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E2010_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E2015_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_BUILT_S_E2020_GLOBE_R2022A_54009_1000_V1_0.tif']
# Select built-up data for Reference year
BTU_REFyear= '\n'.join(s for s in BTU_data_store_info if str(REFyear) in s)
# BTU_REFyear
# Info about built-up data for Reference year
BTU_data_store.describe_data(BTU_REFyear).to_dict()
{'data_id': 'GHS_BUILT_S_E2020_GLOBE_R2022A_54009_1000_V1_0.tif',
'data_type': 'mldataset',
'bbox': [-18041000.0, -9000000.0, 18041000.0, 9000000.0],
'time_range': [None, None],
'dims': {'x': 36082, 'y': 18000},
'spatial_res': 1000.0,
'coords': {'x': {'name': 'x', 'dtype': 'float64', 'dims': ['x']},
'y': {'name': 'y', 'dtype': 'float64', 'dims': ['y']},
'spatial_ref': {'name': 'spatial_ref',
'dtype': 'int64',
'dims': [],
'attrs': {'crs_wkt': 'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]',
'spatial_ref': 'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]',
'GeoTransform': '-18041000.0 1000.0 0.0 9000000.0 0.0 -1000.0'}}},
'data_vars': {'band_1': {'name': 'band_1',
'dtype': 'uint32',
'dims': ['y', 'x'],
'chunks': [512, 512],
'attrs': {'AREA_OR_POINT': 'Area',
'_FillValue': 4294967295,
'scale_factor': 1.0,
'add_offset': 0.0,
'grid_mapping': 'spatial_ref'}}},
'attrs': {'source': 's3://xcube-dcfs/GHSL/GHS-BUILT-S/GHS_BUILT_S_E2020_GLOBE_R2022A_54009_1000_V1_0.tif'},
'num_levels': 6}
# # Info about population data for Reference year
# print(BTU_data_store.open_data(BTU_REFyear).num_levels, '\n')
# print(BTU_data_store.open_data(BTU_REFyear).datasets, '\n')
# BTU_data_store.open_data(BTU_REFyear).grid_mapping
level = 0 # this is the highest resolution, which corresponds to 1km
# Build-up dataset for selected Reference year - clipped by AOI
BTU_ds_0 = BTU_data_store.open_data(BTU_REFyear).get_dataset(level).rio.clip(AOI54009.geometry, AOI54009.crs)
BTU_ds_0 = BTU_ds_0.rename({"band_1": "built-up"})
#BTU_ds_0
# # Info AOI coordinates in the world projection Mollweide (ESRI:54009)
# print(BTU_ds_0.x.min().values, BTU_ds_0.x.max().values)
# print(BTU_ds_0.y.min().values, BTU_ds_0.y.max().values)
POP_data_store = new_data_store("s3",
root="xcube-dcfs/GHSL/GHS-POP/",
storage_options=dict(anon=True)
)
OP_data_store info:
POP_data_store_info=list(POP_data_store.get_data_ids())
POP_data_store_info
['GHS_POP_E1975_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E1980_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E1985_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E1990_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E1995_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E2000_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E2005_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E2010_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E2015_GLOBE_R2022A_54009_1000_V1_0.tif', 'GHS_POP_E2020_GLOBE_R2022A_54009_1000_V1_0.tif']
# Select population data for Reference year
POP_REFyear= '\n'.join(s for s in POP_data_store_info if str(REFyear) in s)
#POP_REFyear
# Info about built-up data for Reference year
POP_data_store.describe_data(POP_REFyear).to_dict()
{'data_id': 'GHS_POP_E2020_GLOBE_R2022A_54009_1000_V1_0.tif',
'data_type': 'mldataset',
'bbox': [-18041000.0, -9000000.0, 18041000.0, 9000000.0],
'time_range': [None, None],
'dims': {'x': 36082, 'y': 18000},
'spatial_res': 1000.0,
'coords': {'x': {'name': 'x', 'dtype': 'float64', 'dims': ['x']},
'y': {'name': 'y', 'dtype': 'float64', 'dims': ['y']},
'spatial_ref': {'name': 'spatial_ref',
'dtype': 'int64',
'dims': [],
'attrs': {'crs_wkt': 'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]',
'spatial_ref': 'PROJCS["World_Mollweide",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mollweide"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]',
'GeoTransform': '-18041000.0 1000.0 0.0 9000000.0 0.0 -1000.0'}}},
'data_vars': {'band_1': {'name': 'band_1',
'dtype': 'float64',
'dims': ['y', 'x'],
'chunks': [512, 512],
'attrs': {'AREA_OR_POINT': 'Area',
'STATISTICS_COVARIANCES': 330985.3282731281,
'STATISTICS_MAXIMUM': 615756.8270821,
'STATISTICS_MEAN': 56.450798856297,
'STATISTICS_MEDIAN': 0.0,
'STATISTICS_MINIMUM': 0,
'STATISTICS_SKIPFACTORX': 1,
'STATISTICS_SKIPFACTORY': 1,
'STATISTICS_STDDEV': 575.31324361006,
'STATISTICS_VALID_PERCENT': 21.26,
'_FillValue': -200.0,
'scale_factor': 1.0,
'add_offset': 0.0,
'grid_mapping': 'spatial_ref'}}},
'attrs': {'source': 's3://xcube-dcfs/GHSL/GHS-POP/GHS_POP_E2020_GLOBE_R2022A_54009_1000_V1_0.tif'},
'num_levels': 6}
# # Info about population data for Reference year
# print(POP_data_store.open_data(POP_REFyear).num_levels, '\n')
# print(POP_data_store.open_data(POP_REFyear).datasets, '\n')
# POP_data_store.open_data(POP_REFyear).grid_mapping
level = 0 # this is the highest resolution, which corresponds to 1km
# Population dataset for selected Reference year - clipped by AOI
POP_ds_0 = POP_data_store.open_data(POP_REFyear).get_dataset(level).rio.clip(AOI54009.geometry, AOI54009.crs)
POP_ds_0 = POP_ds_0.rename({"band_1": "population_grid"})
#POP_ds_0
BTU_da=(BTU_ds_0['built-up']).squeeze()
BTU_da = BTU_da.load()
BTU_da=BTU_da.where(BTU_da != BTU_da.attrs['_FillValue']).fillna(0)
BTU_da
<xarray.DataArray 'built-up' (y: 831, x: 1266)>
array([[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[1203., 0., 0., ..., 0., 0., 0.],
...,
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.]])
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_ref
_FillValue: 4294967295POP_da=(POP_ds_0['population_grid']).squeeze()
POP_da = POP_da.load()
POP_da=POP_da.where(POP_da != POP_da.attrs['_FillValue']).fillna(0)
POP_da
<xarray.DataArray 'population_grid' (y: 831, x: 1266)>
array([[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[5.30303578, 0. , 0. , ..., 0. , 0. ,
0. ],
...,
[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[0. , 0. , 0. , ..., 0. , 0. ,
0. ],
[0. , 0. , 0. , ..., 0. , 0. ,
0. ]])
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0
Attributes:
AREA_OR_POINT: Area
STATISTICS_COVARIANCES: 330985.3282731281
STATISTICS_MAXIMUM: 615756.8270821
STATISTICS_MEAN: 56.450798856297
STATISTICS_MEDIAN: 0.0
STATISTICS_MINIMUM: 0
STATISTICS_SKIPFACTORX: 1
STATISTICS_SKIPFACTORY: 1
STATISTICS_STDDEV: 575.31324361006
STATISTICS_VALID_PERCENT: 21.26
scale_factor: 1.0
add_offset: 0.0
grid_mapping: spatial_ref
_FillValue: -200.0# # Info about DataArray
# print(type(POP_da),'\n')
# print(type(POP_da.data))
# print(POP_da.name), print()
# print("The shape of data is:\n", POP_da.shape,'\n')
# print(POP_da.coords), print()
# print("the minimum raster value is:\n", np.nanmin(POP_da.values))
# print("the maximum raster value is:\n", np.nanmax(POP_da.values),'\n')
# print('dim:',POP_da.dims,'\n')
# print('attributes/metadata:\n',POP_da.attrs,'\n')
# print('----------------')
# # View generate metadata associated with the raster file
# # World Mollweide ESRI:54009,
# print("The crs of data is:\n", POP_da.rio.crs,'\n')
# print('transform:\n',POP_da.rio.transform(),'\n')
# print("The spatial extent is:\n", POP_da.rio.bounds(),'\n')
# print("The nodatavalue of data is:\n", POP_da.rio.nodata,'\n')
# print("The spatial resolution of data is:\n", POP_da.rio.resolution(),'\n')
# print(BTU_da.data.min(), BTU_da.data.max())
# print(POP_da.data.min(), POP_da.data.max().round(0))
# POP_da.data
# HDC
# built-up area
HDC_BTU_mask = 1 * np.ones_like(BTU_da) * (BTU_da > HDC_Bdens).astype(np.int8) # HDC_Bdens (share of built-up area/km2: 0.5 => 500000 m2/km2)
# population
HDC_POP_mask = 1 * np.ones_like(POP_da) * (POP_da > HDC_Pdens).astype(np.int8) # HDC_Pdens >=1500 inh/km2
HDC_UC01_mask = HDC_BTU_mask + HDC_POP_mask
# MDC
MDC_POP_mask = 1 * np.ones_like(POP_da) * (POP_da > MDC_Pdens).astype(np.int8) # MDC_Pdens >=300 inh/km2
# BINARY MASK:
HDC_UC01_mask = 1 * np.ones_like(HDC_UC01_mask) * (HDC_UC01_mask> 0).astype(np.int8)
MDC_UC01_mask = MDC_POP_mask
https://scikit-image.org/docs/dev/api/skimage.measure.html#skimage.measure.label
# 4/8-connectivity regions:
# HDC_Tcon=1 for 4-connectivity clusters
# MDC_Tcon=2 for 8-connectivity clusters
HDC_UC02_4conn = measure.label(HDC_UC01_mask, background=0, return_num=False, connectivity=HDC_Tcon)
# print('HDC number of regions:', HDC_UC02_4conn.max())
MDC_UC02_8conn = measure.label(MDC_UC01_mask, background=0, return_num=False, connectivity=MDC_Tcon)
# print('MDC number of regions:', MDC_UC02_8conn.max())
HDC_UC03=np.zeros_like(HDC_POP_mask).astype('int')
POPregionIDsum=[]
POPregionSUM=[]
for regionID in range(1, HDC_UC02_4conn.max()):
POPregionID=((HDC_UC02_4conn == regionID).astype('int')) * POP_da.data
# threshold for population in cluster > 50000
if POPregionID.sum() > HDC_Pmin:
POPregionSUM.append(int(round(POPregionID.sum())))
POPregionIDsum.append(regionID)
HDC_UC03=HDC_UC03+np.where(HDC_UC02_4conn == regionID, HDC_UC02_4conn, 0)
MDC_UC03=np.zeros_like(MDC_POP_mask).astype('int')
POPregionIDsum=[]
POPregionSUM=[]
for regionID in range(1, MDC_UC02_8conn.max()):
POPregionID=((MDC_UC02_8conn == regionID).astype('int')) * POP_da.data
# threshold for population in cluster > 5000
if POPregionID.sum() > MDC_Pmin:
POPregionSUM.append(int(round(POPregionID.sum())))
POPregionIDsum.append(regionID)
MDC_UC03=MDC_UC03+np.where(MDC_UC02_8conn == regionID, MDC_UC02_8conn, 0)
inverse raster / gaps
HDC_UC03_BIN=np.where(HDC_UC03 > 0, 1, 0)
HDC_UC03_INV=np.where(HDC_UC03_BIN == 1, 0, 1)
HDC_UC03_INV = measure.label(HDC_UC03_INV, background=0, connectivity=2) # 8-connectivity
# number of holes
# HDC_UC03_INV.max()
gaps filling
unique, counts = np.unique(HDC_UC03_INV, return_counts=True,)
gaps=list(zip(unique, counts))
#print(gaps)
for gap in gaps:
if gap[1] < HDC_GAPfill:
gapID=(HDC_UC03_INV == gap[0]).astype(int)
HDC_UC03_BIN=HDC_UC03_BIN+gapID
regions where POP > 50000 with filled gaps < 15km2
HDC_UC03_4conn = measure.label(HDC_UC03_BIN, background=0, connectivity=HDC_Tcon) # connectivity=HDC_Tcon
print(HDC_UC03_4conn.min(),HDC_UC03_4conn.max())
0 123
# binary raster
HDC_UC04_BIN=np.where(HDC_UC03_4conn >= 1, 1, 0)
print(HDC_UC04_BIN.sum())
8165
closing, erosion, dilatation
HDC_UC04=np.zeros_like(HDC_UC03_BIN).astype('int')
# regiony sa vyhodnocuju sekvencne v cykle:
POPregionIDsum=[]
POPregionSUM=[]
for regionID in range(1, HDC_UC03_4conn.max()):
# mask 0-1 for region, # POP in region
POPregionID=((HDC_UC03_4conn == regionID).astype('int'))
#print(POPregionID.min(), POPregionID.max(), POPregionID.sum())
POPregionID = morphology.binary_closing(POPregionID, square(3))
POPregionID= morphology.binary_dilation(POPregionID)
POPregionID= morphology.binary_erosion(POPregionID)
#POPregionID = skimage.morphology.binary_closing(POPregionID, square(3))
HDC_UC04=HDC_UC04+POPregionID
#print(HDC_UC05.sum())
# gaps filling after morphology:
HDC_UC04_INV=np.where(HDC_UC04 == 1, 0, 1)
HDC_UC04_INV = measure.label(HDC_UC04_INV, background=0, connectivity=2)
# number of holes
HDC_UC04_INV.max()
2
unique1, counts1 = np.unique(HDC_UC04_INV, return_counts=True,)
gaps1 = [x for x in zip(unique1, counts1) if x[1] < HDC_GAPfill]
#print(gaps1)
for gap1 in gaps1:
gapID1=(HDC_UC04_INV == gap1[0]).astype(int)
HDC_UC04=HDC_UC04+gapID1
HDC_UC04_4conn = measure.label(HDC_UC04, background=0, connectivity=HDC_Tcon)
HDC_UC04_4conn.max()
121
# binary raster
HDC_UC04_4conn_BIN=np.where(HDC_UC04_4conn >= 1, 1, 0)
HDC_UC04_4conn_BIN.sum()
9193
# binary
MDC_UC03=np.where(MDC_UC03 > 0, 1, 0)
#
MDC_UC04=MDC_UC03-HDC_UC04
MDC_UC04=np.where(MDC_UC04 == 1, 1, 0)
MDC_UC04_8conn = measure.label(MDC_UC04, background=0, connectivity=MDC_Tcon)
print(MDC_UC04_8conn.min(),MDC_UC04_8conn.max())
0 1909
HDC and MDC to vector
HDC: numpy.ndarray to DataArray
HDC_UC05 = xr.DataArray(
data=HDC_UC04_4conn,
dims=BTU_da.dims,
coords=BTU_da.coords,
name='urban_centre',
)
HDC_UC05=HDC_UC05.astype(np.int16)
HDC_UC05
<xarray.DataArray 'urban_centre' (y: 831, x: 1266)>
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=int16)
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0# # Info about HDC DataArray
# print('type:',type(HDC_UC05))
# print('name:',HDC_UC05.name)
# print('shape:',HDC_UC05.shape)
# #print(HDC_UC05.coords)
# print('dim:',HDC_UC05.dims,'\n')
# # World Mollweide ESRI:54009
# print("The crs is: ", HDC_UC05.rio.crs)
# print('transform: ',HDC_UC05.rio.transform())
# print("The spatial extent is: ", HDC_UC05.rio.bounds())
# print("The nodatavalue of data is: ", HDC_UC05.rio.nodata)
# print("The shape of data is: ", HDC_UC05.shape)
# print("The spatial resolution of data is: ", HDC_UC05.rio.resolution())
# print("the minimum raster value is: ", np.nanmin(HDC_UC05.values))
# print("the maximum raster value is: ", np.nanmax(HDC_UC05.values))
# print("The metadata of data is: ", HDC_UC05.attrs)
HDC: export to vector
# mask=None
mask=HDC_UC05.data != 0
with HDC_UC05 as src:
image = src.data # first band
results = (
{'properties': {'raster_val': v}, 'geometry': s}
for i, (s, v)
in enumerate(
shapes(image, mask=mask, connectivity=4, transform=HDC_UC05.rio.transform())))
geoms = list(results)
# first feature
# print(geoms[0])
HDC_UrbanCentre = gpd.GeoDataFrame.from_features(geoms)
HDC_UrbanCentre.crs=HDC_UC05.rio.crs
HDC_UrbanCentre.raster_val=HDC_UrbanCentre.raster_val.astype('int')
HDC_UrbanCentre.rename(columns={'raster_val':'hdcid'}, inplace=True)
HDC_UrbanCentre['area_km2'] = HDC_UrbanCentre['geometry'].area/ 10**6
HDC_UrbanCentre['code']=3
HDC_UrbanCentre.sort_values(['area_km2'],ascending=False).head(3)
| geometry | hdcid | area_km2 | code | |
|---|---|---|---|---|
| 32 | POLYGON ((991000.000 6151000.000, 994000.000 6... | 29 | 729.0 | 3 |
| 70 | POLYGON ((1423000.000 5916000.000, 1424000.000... | 69 | 478.0 | 3 |
| 39 | POLYGON ((1542000.000 6118000.000, 1551000.000... | 36 | 457.0 | 3 |
MDC: numpy.ndarray to DataArray
MDC_UC05 = xr.DataArray(
data=MDC_UC04_8conn,
dims=BTU_da.dims,
coords=BTU_da.coords,
name='urban_cluster',
)
MDC_UC05=MDC_UC05.astype(np.int16)
MDC_UC05
<xarray.DataArray 'urban_cluster' (y: 831, x: 1266)>
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=int16)
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0# # BASIC INFO
# print('type:',type(MDC_UC05))
# print('name:',MDC_UC05.name)
# print('shape:',MDC_UC05.shape)
# #print(MDC_UC05.coords)
# print('dim:',MDC_UC05.dims,'\n')
# # View generate metadata associated with the raster file
# # World Mollweide ESRI:54009
# print("The crs is: ", MDC_UC05.rio.crs)
# print('transform: ',MDC_UC05.rio.transform())
# print("The spatial extent is: ", MDC_UC05.rio.bounds())
# print("The nodatavalue of data is: ", MDC_UC05.rio.nodata)
# print("The shape of data is: ", MDC_UC05.shape)
# print("The spatial resolution of data is: ", MDC_UC05.rio.resolution())
# print("the minimum raster value is: ", np.nanmin(MDC_UC05.values))
# print("the maximum raster value is: ", np.nanmax(MDC_UC05.values))
# print("The metadata of data is: ", MDC_UC05.attrs,'\n')
# print(MDC_UC05.data.max())
MDC: export to vector
# mask=None
mask=MDC_UC05.data != 0
with MDC_UC05 as src:
image = src.data # first band
results = (
{'properties': {'raster_val': v}, 'geometry': s}
for i, (s, v)
in enumerate(
shapes(image, mask=mask, connectivity=8, transform=MDC_UC05.rio.transform())))
geoms = list(results)
MDC_UrbanCluster = gpd.GeoDataFrame.from_features(geoms)
MDC_UrbanCluster.crs=MDC_UC05.rio.crs
MDC_UrbanCluster.raster_val=MDC_UrbanCluster.raster_val.astype('int')
MDC_UrbanCluster.rename(columns={'raster_val':'mdcid'}, inplace=True)
MDC_UrbanCluster["area_km2"] = MDC_UrbanCluster['geometry'].area/ 10**6
MDC_UrbanCluster['code']=2
MDC_UrbanCluster.sort_values(['area_km2'],ascending=False).head(3)
| geometry | mdcid | area_km2 | code | |
|---|---|---|---|---|
| 1164 | POLYGON ((1407000.000 5892000.000, 1408000.000... | 997 | 849.0 | 2 |
| 505 | POLYGON ((1568000.000 6122000.000, 1569000.000... | 423 | 632.0 | 2 |
| 1141 | POLYGON ((1437000.000 5876000.000, 1438000.000... | 1047 | 538.0 | 2 |
# Create folder for GHSL_results, if it doesn't exist:
Path(os.path.realpath('GHSL_results')).mkdir(parents=True, exist_ok=True)
# EXPORT 2 shapefile World Mollweide (ESRI:54009)
HDC_UrbanCentre.to_file(os.path.join(os.path.realpath('GHSL_results'),'HDC_UrbanCentre_'+str(REFyear)+'.shp'))
# Export to shapefile WGS84 (EPSG:4326 )
# HDC_UrbanCentre.to_crs(epsg=4326).to_file(os.path.join(os.path.realpath('GHSL_results'),'HDC_UrbanCentre_4326_'+str(REFyear)+'.shp'))
# EXPORT 2 shapefile World Mollweide (ESRI:54009)
MDC_UrbanCluster.to_file(os.path.join(os.path.realpath('GHSL_results'),'MDC_UrbanCluster_'+str(REFyear)+'.shp'))
# Export to shapefile WGS84 (EPSG:4326 )
# MDC_UrbanCentre.to_crs(epsg=4326).to_file(os.path.join(os.path.realpath('GHSL_results'),'MDC_UrbanCluster_4326_'+str(REFyear)+'.shp'))
fig = Figure(width=600, height=600)
m = AOI54009.explore(style_kwds={'fillColor':'None','color':'blue','weight': 1}, name='aoi')
# parallels:
for lat in range(-90, 91, lat_interval):
folium.PolyLine([[lat, -180],[lat, 180]], weight=0.5).add_to(m)
# meridianss:
for lon in range(-180, 181, lon_interval):
folium.PolyLine([[-90, lon],[90, lon]], weight=0.5).add_to(m)
MDC_UrbanCluster.explore(m=m,style_kwds={'fillColor':'orange','color':'orange','weight': 1}, name='MDC_UrbanCluster')
HDC_UrbanCentre.explore(m=m,style_kwds={'fillColor':'red','color':'red','weight': 1.5}, name='HDC_UrbanCentre')
folium.LayerControl('topright',collapsed=True).add_to(m)
minimap = plugins.MiniMap(height=90, width=90, zoom_level_offset=-7)
m.add_child(minimap)
fig.add_child(m)
m
# MASK
DUC_mask=HDC_UC01_mask
#print(DUC_mask.sum())
# CONNECTIVITY
DUC_4conn = measure.label(DUC_mask, background=0, return_num=False, connectivity=HDC_Tcon) #.squeeze()
print('DUC number of regions:', DUC_4conn.max())
DUC number of regions: 2437
# TOTAL POPULATION IN CLUSTER
DUC_UC03=np.zeros_like(HDC_POP_mask).astype('int')
POPregionIDsum=[]
POPregionSUM=[]
for regionID in range(1, DUC_4conn.max()):
POPregionID=((DUC_4conn == regionID).astype('int')) * POP_da.data
# threshold for population in cluster > 5000
if POPregionID.sum() > MDC_Pmin:
POPregionSUM.append(int(round(POPregionID.sum())))
POPregionIDsum.append(regionID)
DUC_UC03=DUC_UC03+np.where(DUC_4conn == regionID, DUC_4conn, 0)
# SUBTRACTION OF HDC
DUC_UC03_BIN=np.where(DUC_UC03 > 0, 1, 0)
DUC_UC03_BIN=DUC_UC03_BIN-HDC_UC04_4conn_BIN
DUC_UC03_BIN=np.where(DUC_UC03_BIN >0, 1, 0)
print(DUC_UC03_BIN.sum())
4591
# CONNECTIVITY
DUC_UC04_4conn = measure.label(DUC_UC03_BIN, background=0, connectivity=HDC_Tcon) # connectivity=HDC_Tcon
print(DUC_UC04_4conn.min(),DUC_UC04_4conn.max())
0 756
# BINARY MASK
DUC_UC04_4conn_BIN=np.where(DUC_UC04_4conn > 0, 1, 0)
DUC_UC04_4conn_BIN.sum()
4591
DUC: numpy.ndarray to DataArray
DUC_UC05 = xr.DataArray(
data=DUC_UC04_4conn,
dims=BTU_da.dims,
coords=BTU_da.coords,
name='dense_urban_cluster',
)
DUC_UC05=DUC_UC05.astype(np.int16)
DUC_UC05
<xarray.DataArray 'dense_urban_cluster' (y: 831, x: 1266)>
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=int16)
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0DUC: export to vector
# mask=None
mask=DUC_UC05.data != 0
with DUC_UC05 as src:
image = src.data # first band
results = (
{'properties': {'raster_val': v}, 'geometry': s}
for i, (s, v)
in enumerate(
shapes(image, mask=mask, connectivity=4, transform=DUC_UC05.rio.transform())))
geoms = list(results)
DUC_DenseUrbanCluster = gpd.GeoDataFrame.from_features(geoms)
DUC_DenseUrbanCluster.crs=DUC_UC05.rio.crs
DUC_DenseUrbanCluster.raster_val=DUC_DenseUrbanCluster.raster_val.astype('int')
DUC_DenseUrbanCluster.rename(columns={'raster_val':'ducid'}, inplace=True)
DUC_DenseUrbanCluster['code']=23
# MASK
SDUC_mask=MDC_UC01_mask
print(DUC_mask.sum())
<xarray.DataArray ()>
array(14494.)
Coordinates:
spatial_ref int64 0
# CONNECTIVITY
SDUC_8conn = measure.label(SDUC_mask, background=0, return_num=False, connectivity=MDC_Tcon) # 8-connectivity cluster
print('SDUC number of regions:', SDUC_8conn.max())
SDUC number of regions: 15569
# TOTAL POPULATION IN CLUSTER
SDUC_UC03=np.zeros_like(HDC_POP_mask).astype('int')
POPregionIDsum=[]
POPregionSUM=[]
for regionID in range(1, SDUC_8conn.max()):
POPregionID=((SDUC_8conn == regionID).astype('int')) * POP_da.data
# threshold for population in cluster > 5000
if POPregionID.sum() > MDC_Pmin:
POPregionSUM.append(int(round(POPregionID.sum())))
POPregionIDsum.append(regionID)
SDUC_UC03=SDUC_UC03+np.where(SDUC_8conn == regionID, SDUC_8conn, 0)
# Minus HDC and DUC
SDUC_UC03_BIN=np.where(SDUC_UC03 > 0, 1, 0)
SDUC_UC03_BIN=SDUC_UC03_BIN - HDC_UC04_4conn_BIN - DUC_UC04_4conn_BIN
SDUC_UC03_BIN=np.where(SDUC_UC03_BIN >0, 1, 0)
print(SDUC_UC03_BIN.sum())
32703
SDUC_UC04_8conn = measure.label(SDUC_UC03_BIN, background=0, connectivity=MDC_Tcon) # connectivity=HDC_Tcon
print(SDUC_UC04_8conn.min(),SDUC_UC04_8conn.max())
0 2536
SDUC_UC05 = xr.DataArray(
data=SDUC_UC04_8conn,
dims=BTU_da.dims,
coords=BTU_da.coords,
name='semidense_urban_cluster',
)
SDUC_UC05=SDUC_UC05.astype(np.int16)
SDUC_UC05
<xarray.DataArray 'semidense_urban_cluster' (y: 831, x: 1266)>
array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=int16)
Coordinates:
* x (x) float64 7.085e+05 7.095e+05 ... 1.972e+06 1.974e+06
* y (y) float64 6.386e+06 6.386e+06 ... 5.558e+06 5.556e+06
spatial_ref int64 0SDUC: export to vector
# mask=None
mask=SDUC_UC05.data != 0
with SDUC_UC05 as src:
image = src.data # first band
results = (
{'properties': {'raster_val': v}, 'geometry': s}
for i, (s, v)
in enumerate(
shapes(image, mask=mask, connectivity=8, transform=SDUC_UC05.rio.transform())))
geoms = list(results)
# first feature
# print(geoms[0])
SDUC_SemiDenseUrbanCluster_temp = gpd.GeoDataFrame.from_features(geoms)
SDUC_SemiDenseUrbanCluster_temp.crs=SDUC_UC05.rio.crs
SDUC_SemiDenseUrbanCluster_temp.raster_val=SDUC_SemiDenseUrbanCluster_temp.raster_val.astype('int')
SDUC_SemiDenseUrbanCluster_temp.rename(columns={'raster_val':'sducid'}, inplace=True)
SDUC_SemiDenseUrbanCluster_temp['code']=22
# HDC & DUC - BUFFER 3km
res_union_buffer = gpd.overlay(HDC_UrbanCentre[['geometry']], DUC_DenseUrbanCluster[['geometry']], how='union',keep_geom_type=True).buffer(3000).unary_union
res_union_buffer =gpd.GeoSeries(res_union_buffer)
res_union_buffer = gpd.GeoDataFrame(geometry=gpd.GeoSeries(res_union_buffer))
res_union_buffer.crs=HDC_UrbanCentre.crs
res_union_buffer[['buff3']]=1
# sjoin with buffer
SDUC_SemiDenseUrbanCluster_temp = gpd.sjoin(SDUC_SemiDenseUrbanCluster_temp, res_union_buffer[['buff3', 'geometry']], how='left',predicate='intersects')
SDUC_SemiDenseUrbanCluster_temp.drop(columns =['index_right'],inplace=True)
SDUC_SemiDenseUrbanCluster_temp['code'] = np.where(SDUC_SemiDenseUrbanCluster_temp['buff3'] == 1, 21, 22)
# print(SDUC_SemiDenseUrbanCluster_temp.tail(2))
# SDUC_SemiDenseUrbanCluster
SDUC_SemiDenseUrbanCluster = SDUC_SemiDenseUrbanCluster_temp[SDUC_SemiDenseUrbanCluster_temp.code == 22].drop(columns=['buff3'])
# PUC_PeriUrbanCluster
PUC_PeriUrbanCluster = SDUC_SemiDenseUrbanCluster_temp[SDUC_SemiDenseUrbanCluster_temp.code == 21].drop(columns=['buff3'])
PUC_PeriUrbanCluster.rename(columns={'sducid':'pucid'},inplace=True)
DUC: export to shapefile
# EXPORT 2 shapefile World Mollweide (ESRI:54009)
DUC_DenseUrbanCluster.to_file(os.path.join(os.path.realpath('GHSL_results'),'DUC_DenseUrbanCluster_'+str(REFyear)+'.shp'))
# Export to shapefile WGS84 (EPSG:4326 )
# DUC_DenseUrbanCluster.to_crs(epsg=4326).to_file(os.path.join(os.path.realpath('GHSL_results'),'DUC_DenseUrbanCluster_4326_'+str(REFyear)+'.shp'))
SDUC: export to shapefile
# EXPORT 2 shapefile World Mollweide (ESRI:54009)
SDUC_SemiDenseUrbanCluster.to_file(os.path.join(os.path.realpath('GHSL_results'),'SDUC_SemiDenseUrbanCluster_'+str(REFyear)+'.shp'))
# Export to shapefile WGS84 (EPSG:4326 )
# SDUC_SemiDenseUrbanCluster.to_crs(epsg=4326).to_file(os.path.join(os.path.realpath('GHSL_results'),'SDUC_SemiDenseUrbanCluster_4326_'+str(REFyear)+'.shp'))
PUC: export to shapefile
# EXPORT 2 shapefile World Mollweide (ESRI:54009)
PUC_PeriUrbanCluster.to_file(os.path.join(os.path.realpath('GHSL_results'),'PUC_PeriUrbanCluster_'+str(REFyear)+'.shp'))
# Export to shapefile WGS84 (EPSG:4326 )
# PUC_PeriUrbanCluster.to_crs(epsg=4326).to_file(os.path.join(os.path.realpath('GHSL_results'),'PUC_PeriUrbanCluster_4326'+str(REFyear)+'.shp'))
fig = Figure(width=600, height=600)
m = AOI54009.explore(style_kwds={'fillColor':'None','color':'blue','weight': 1}, name='aoi')
# parallels:
for lat in range(-90, 91, lat_interval):
folium.PolyLine([[lat, -180],[lat, 180]], weight=0.5).add_to(m)
# meridianss:
for lon in range(-180, 181, lon_interval):
folium.PolyLine([[-90, lon],[90, lon]], weight=0.5).add_to(m)
PUC_PeriUrbanCluster.explore(m=m,style_kwds={'fillColor':'yellow','color':'yellow','weight': 1.5}, name='PUC_PeriUrbanCluster')
SDUC_SemiDenseUrbanCluster.explore(m=m,style_kwds={'fillColor':'darkgoldenrod','color':'darkgoldenrod','weight': 1.5}, name='SDUC_SemiDenseUrbanCluster')
DUC_DenseUrbanCluster.explore(m=m,style_kwds={'fillColor':'brown','color':'brown','weight': 1.5}, name='DUC_DenseUrbanCluster')
MDC_UrbanCluster.explore(m=m,style_kwds={'fillColor':'orange','color':'orange','weight': 1}, name='MDC_UrbanCluster', show = False)
HDC_UrbanCentre.explore(m=m,style_kwds={'fillColor':'red','color':'red','weight': 1.5}, name='HDC_UrbanCentre')
folium.LayerControl('topright',collapsed=True).add_to(m)
minimap = plugins.MiniMap(height=90, width=90, zoom_level_offset=-7)
m.add_child(minimap)
fig.add_child(m)
m